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Abstract. In this paper, we introduce the convergence analysis of the fixed pivot technique given 
by S. Kumar and Ramkrishna 1251 for the nonlinear aggregation population balance equations which 
are of substantial interest in many areas of science: colloid chemistry, aerosol physics, astrophysics, 
polymer science, oil recovery dynamics, and mathematical biology. In particular, we investigate the 
convergence for five different types of uniform and non-uniform meshes which turns out that the fixed 
pivot technique is second order convergent on a uniform and non-uniform smooth meshes. Moreover, 
it yields first order convergence on a locally uniform mesh. Finally, the analysis exhibits that the 
method does not converge on an oscillatory and non-uniform random meshes. Mathematical results 
of the convergence analysis are also demonstrated numerically. 
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1 . Introduction. The continuous aggregation population balance equation (PBE) 
or Smoluchowski coagulation equation describes the kinetics of particle growth in 
which particles can aggregate via binary interaction to form larger particles. This 
model arises in many fields of science and engineering: kinetics of phase transforma- 
tions in binary alloys such as segregation of binary alloys, aggregation of red blood 
cells in biology, fluidized bed granulation processes, aerosol physics, i.e the evolution 
of a system of solid or liquid particles suspended in a gas, formation of planets in 
astrophysics, polymer science and many more. The nonlinear continuous aggregation 
population balance equation is given by 

(L1) If* R{x _ ^ x _ y)f{ ^ y)dy 



dt 2 



K(x,y)f(t,x)f(t,y)dy, 



with f(x,0) = f m (x) > 0, xe]Q,oo[. 

where the variables x > and t > denote the size of the particles and time respec- 
tively. The number density of particles of size x at time t is denoted by f{x,t) > 0. 
The aggregation kernel K (x, y) > represents the rate at which particles of size x 
coalesce with those of size y. It will be assumed throughout that K{x, y) — K(y, x) 
for all x, y > 0, i.e. symmetric and K (x, y) = for either x = or y = 0. 
Mathematical results on existence and uniqueness of solutions to equation can be 
found in[7l[9l[l2l[T9l[23[Ml[29l[3a[33[38]for different classes of aggregation kernels. 
The pure aggregation PBE (jl.ip can be solved analytically only for some specific 
examples of aggregation kernels, see [8j [TBI KZ]- I" general we need to solve them 
numerically. To apply a numerical method, first we need to consider the following 
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truncated form of the problem Ijl.ip by taking a finite computational domain. 



where n(t,x) represents the solution to the truncated equation (|1.2[) . The existence 
and uniqueness of non-negative solutions for the truncated PBE (jl.2l) has been shown 
in [51 O I3H]- In HH HH1 I3S], it is proven that the sequence of solutions to 
the truncated problems converge weakly to the solution of the original problem in a 
weighted L 1 space as xmax — > oo for certain classes of kernels. 

Many numerical methods have been proposed to solve the truncated aggregation 
PBE CL2): finite element methods p31 [35j [36], finite volume methods [3] [T4] ITS], 
stochastic methods [TO] HU [30] , moment methods [39] and sectional methods [22[ [23] 
I25[ 126] . By implementing most of these methods, we may have a quite satisfactory 
results for the number density but not for moments. However, the moment methods 
give opposite results. To have a satisfactory information for the number density dis- 
tribution as well as some selected moments, the sectional methods have become more 
useful nowadays. Several authors have proposed sectional methods for aggregation 
PBE: S. Kumar and Ramkrishna (25j [SB], J. Kumar et al. [22] [23] and Vanni [40] . 
The fixed pivot technique given by S. Kumar and Ramkrishna is the most extensively 
used sectional method. This technique also efficiently works for a multi-dimensional 
size variable [3]. 

Recently J. Kumar and Warnecke [24] have published the numerical analysis of 
the fixed pivot technique for breakage PBEs. This case was simpler due to the linearity 
of that equation. However the convergence analysis of the technique was still open 
for aggregation PBEs (|1.2[) . This was a challenging task due to the non-linearity of 
the equation. So the purpose of this work is to demonstrate the missing convergence 
analysis of the fixed pivot technique for aggregation PBEs in the literature. 

Let us now briefly outline the contents of this paper. Along with the general idea 
of sectional methods, a concise review of the mathematical formulation of the fixed 
pivot technique is given in the following section. A theorem from Hundsdorfer and 
Verwer |20] used in further analysis and the main result for the convergence of the 
fixed pivot technique are also stated in Section [2] To show the convergence of the 
scheme for solving aggregation PBEs (|1.2[) . the consistency and Lipschitz conditions 
are discussed in Sections [3] and HI respectively. Numerical simulations are performed 
in Section [5J And, Section [5] states some conclusions. 

2. The Sectional Methods. The mathematical formulation of sectional meth- 
ods is reviewed from [53] Section 2]. These methods calculate the total number of 
particles in finite number of cells. First of all, the continuous interval Q :=]0, xrnax] 
is divided into a finite number of cells defining size classes A4 :=]Xi_i/2, £1+1/2], * = 
1,...,/. Set 

X1/2 = 0, x I+ x/ 2 = imax, Ax min - = ^+1/2 - < Az. 

For the purpose of later analysis we assume quasi uniformity of the grids, i.e. 



(1.2) 




with 



n(x, 0) 



n m (x) > 0, x 



£ O :=]0,xmax], 



(2.1) 



Ax 



<C, 



Ax. 
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where C is a positive constant. The center point of each cell Xi = (x i _ 1 ^ 2 + 
x i+1 / 2 )/2, i — 1, ... ,1 is called pivot or grid point. The integration of the truncated 
PBE (|1.2[) over each cell gives a semi- discrete system in M. 1 

dN 

(2.2) = B D, with N(0) = N m . 

at 

where N m , N, B, D £ R 1 . The ith components of vectors N,N m ,B, and D are 
respectively, defined as 

(2.3) Ni(t)= [ 1+1/2 n(t,x)dx, with = [ n in (x)dx, 



(2.4) Bi = - / K(x - y,y)n(t,x - y)n(t,y)dydx. 

1 Jx i _ 1/2 Jo 



and 



r x i+i/2 r x i+i/2 

(2.5) Di= / K(x,y)n(t,y)n(t,x)dydx. 

J' x i-l/2 JO 

Here the vector N is formed by the vector of values of the step function obtained by 
L 2 projection of the exact solution n into the space of step functions constant on each 
cell. Note that this projection error can easily be shown of second order, see section 
5.2.2 in |18j . Finally, by taking numerical approximations of Bi and Di in terms of 
Ni(t), the sectional methods give the following discretized form 

r/N ...... 

(2.6) = B(N) — D(N) =: F(t, N), with N(0) = N m , 

at 

where N, B, D £ M. 1 . The total number of particles in ith cell, Ni(t) are numeri- 
cally approximated by the ith component, Ni(i) of the vector N. The ith components 
of B and D arc denoted by Bi and Di respectively which are defined in spatially dis- 
cretized system (|2.9[) obtained by the fixed pivot technique. 

2.1. The Fixed Pivot Technique. The fixed pivot technique relies on the 
following idea of birth modification. In , the truncated PBE (|1.2p is modified to 

(2.7) — / n(t,x)dx^- A+(x) / K(x - y,y)n(t,x - y)n(t,y)dydx 



dt L 2 

■ -1/2 







K( x ) / K(x - y,y)n(t,x - y)n(t,y)dydx 
x i-i Jo 
x i+i/2 r x i+i/2 

/ K(x,y)n(t,y)n(t,x)dydx. 

x i-l/2 JO 

where 

(2.8) Xf(x)= x ~ Xi ±\ i = l,...,I. 

Xi - Xi±l 

According to S. Kumar and Ramkrishna in |25j . the first and second integral terms 
on the right hand side in (I2.7[) set to be zero for i = / and i = 1, respectively. 
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Inserting the number density approximation n(t,x) ss Y]j—i Ni(t)6(x — £Ej), into 
the above equation, we get the following spatially discretized system 



(2.9) —jT= E U-^Ax+ixk + x^Kix^x^N^Nkit) 

j>k / 1 \ 

E 1 - 9^< fc A ^ + x^Kix^x^Nj^Nkit) 



Xi <Xj-\-Xk <Xi + i 

+ E 

=Bi — Di, 

where Bi and Di represent the discretized birth and death terms, respectively in the 
«th cell obtained from the fixed pivot technique. The basic idea of the fixed pivot 
technique can be described as follows. Assume that a new born particle of a size, 
which is not positioned at a pivot point of any cell, appears due to the aggregation of 
two smaller particles. The particle has to be assigned onto neighboring pivot points 
in such a way that the particle number and mass are conserved. This problem can 
be solved in a unique way. The resulting technique gives very often quite satisfactory 
results. However, the undesirable part is that the fixed pivot technique turns into 
a zero order method on oscillatory and non-uniform random meshes for aggregation 
problems. 

It is important to mention here that at the last boundary cell S. Kumar and 
Ramkrishna simply set the first integral on the right hand side in (|2.7p to be zero 
in their numerical computations. However, we have observed in our analysis for the 
aggregation problem that this setting at the end boundary cell reduces by one order 
the accuracy of the fixed pivot technique on uniform and non-uniform smooth meshes. 
To overcome this problem, we take an extra grid point at a Axi distance away 
from the grid point x\. In the computations contributions which are larger than xi are 
distributed to xi and aJj+i- This idea is used in Lemma 13. II A similar modification 
should also be used at the first boundary cell for breakage problems in Lemma 3.2 
of [23]. The present form of Lemma 3.2 in [23] is not correct at the first boundary 
cell and also reduces by one order accuracy of the fixed pivot technique for breakage 
PBEs on uniform and non-uniform smooth meshes. 

It should be pointed out here that in this work we consider the following discrete 
L\ norm 1 

l|N||=X>il- 

i=l 

We consider C 2 ([a, b]) as a space of twice continuously diffcrentiablc functions on ]a, b[. 
Note that for the sake of simplicity in our analysis we assume that the aggregation 
kernel satisfies 

(2.10) ^eC 2 (]0,a; ma x]x]0,x m ax]). 

As we move on to the subsequent sections, it will be helpful to revisit some 
definitions and an existing theorem given in Hundsdorfer and Verwer [20] . These will 
be of use while discussing in detail the consistency and the convergence of the fixed 
pivot technique. 
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Let || • || denote any norm on R 7 . 

Definition 2.1. The spatial truncation error is defined by the residual left 
by substituting the exact solution N(i) into equation \2. 61) as 

(2.11) a(t) = ^ - (B (N(t)) - D (N(t)) 



The scheme 112. 6\) is called consistent of order p if, for Ax — > 0, 

\\a(t)\\ = 0(Ax p ), uniformly for allt,Q<t< T. 

Definition 2.2. The global discretization error is defined by 
(2.12) e(t) = N(t) - N(t). 

The scheme \2. 6\) is called convergent of order p if for Ax — > 0, 

||e(t)|| = 0(Ax p ), uniformly for all t, < t < T. 

It is important that the solution obtained by the fixed pivot technique remains 
non-negative for all times. This can be easily shown by using the next well known 
theorem. In the following theorem we write M > for a vector M <E R 7 if all of its 
components are non-negative. 

Theorem 2.3. (Hundsdorfer and Verwer JMI)- Suppose that F(i,M) defined in 
\2. 6)) is continuous and satisfies the Lipschitz condition as 

||F(t,P)-F(i,M)|| <£||P-M|| for all P,MeM ; . 

Then the solution of the semi-discrete system \2. 6\) is non-negative if and only if for 
any vector Mel' with M > we have for any i = 1, . . . , / and allt>0 that Mi — 
implies F t (t,M) > 0. 

Proof. The proof can be found in [20], Chap. 1, Theorem 7.1. □ 

Now we shall state the main result which helps us to show the convergence of the 
fixed pivot technique. 

Theorem 2.4. Let us assume that the Lipschitz conditions on B(N(t)) and 
D(N(t)) are satisfied for < t < T and for all N, N € l' where N and N are the 
projected exact and numerical solutions defined in i2.2\) and \2. 6\) respectively. Then 
a consistent discretization method is also convergent and the convergence is of the 
same order as the consistency. 

Proof. Using the equations (|2~TT1) and (j2~T2")l we have for e(t) = N(t) - N(t) 

|e(t) = a(t) + (B(N) - B(N)) - (D(N) - D(N)). 

We then take the norm on both sides to get 

j t \\e(t)\\ < \\a(t)\\ + ||(B(N) - B(N))|| + ||(D(N) - D(N))||. 

Integrating with respect to t with e(0) = and using the Lipschitz conditions 
we obtain the estimates 

||e(t)|| < [ t \\o-(T)\\dT + 2L f \\e(r)\\dr, 
Jo Jo 
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From this it follows by Gronwall's lemma that 



(2.13) 
where 



|e(i)||<||[e X p(2L*)-l], 



eh = max ||cr(t)| 

0<t<T 



If the scheme is consistent then lim/,^0 &h = 0. This completes the proof of the 
Theorem [231 □ 

Note that the proof of the Theorem 12.41 is motivated by a convergence result in 
Linz [31]. 

To fulfill the requirements of Theorem 12.41 for the convergence of the fixed pivot 
technique we need to prove that the scheme is consistent as well as the birth B(N(i)) 
and death D(N(i)) terms satisfy the Lipschitz conditions. 

3. Consistency. We need the following lemma to investigate the consistency of 
the fixed pivot technique for aggregation PBEs (11.21) which is the same as Lemma 3.2 
in [24] except at the boundary cells. 

Lemma 3.1. Consider a function f £ C 2 ([0, xmax}) and a cell centered parti- 
tioning of the domain [0,xmax] given as — X\—\i% < ... < < 2^+1/2 < 
... < = xmax with pivot points Xi — (2^-1/2 + x i+1 / 2 )/2 and a bound 
Ax > Axi — (x i+ i/2 — Xi-1/2) for all i. IfX^(x) and (x) are given by the definition 
then the following expressions can be obtained for the modification error 



Uf) 



f(x)dx 



-1/2 

f{Xi) 



\+{x)f(x)dx 



\ (x)f(x)dx 



2 

f'{Xi) 



Axi 



( Axi-x + Ax l+ i 



12 



V 2 

(Aa; 4+ i - Axi~i)\ Ax 



( Axj-x + Ax i+ i 
V 2 



+C(Aa; 3 ), for i = 2,..., I- 1, 



Xi + l 



f(x)dx- Xt(x)f{x)dx- \r(x)f{x)dx 

^i — 1/2 ** x i J Xi—i 

^[Axt-Ax^+OiAx 2 ), for i = I, 



r x i+i/2 /-Xi+i 
Mf)= / f(x)dx- / \+(x)f{x)dx 

Jx i _ 1/2 J Xi 

= £^-[3Axi- Ax l+1 ]+0{Ax 2 ), for i = l. 

Note that is the extra grid point introduced in section f2. II 

Proof. First, we consider the modification error 3i for i = 2, . . . I — 1 as follows 

3i(f)= f 1+1/2 f(x)dx- ( f^ +1 X+(x)f(x)dx+ p \r(x)f(x)dx 
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Taylor series expansion of f(x) about Xj in J$ yields 



Ax,. - 



\~l(x)dx+ I \ (x)dx 



Xi+i 



Mf) =f(xi) 

-fix 

Substituting the values of A^ and A~ from (|2.8p into the preceding equation, we obtain 



At (a;) (a; - x l )dx + X t (x)(x - x,)dx ) + 0{Ax 3 ) 



(3.1) 3i(f) =f(xi)[Axi - -(x i+ i - afi-i)] 



(x i+ i - Xi-l){(x i+ i - Xi) - (Xi - Xi-i)} 



0(Ax 3 ). 



For the cell centered grids, i.e. xi = (xj_i/ 2 + x i+1 / 2 )/2, the equation (|3.1I) becomes 



Mf) =^-^-\Ax l - ~(Ax i+ i + Axi-i)] 
12 



{Axi+i - Axi-i)< Axi + 



Axi_ x + Ax i+1 



0{Ax 3 ). 



The expression of modification error at last boundary cell i.e. i = I is given by 



Mf) 



f(x)da 



xi+i 



Xf (x)f(x)dx — / X I (x)f(x)dx. 



Applying the midpoint, left rectangle and right rectangle rules in first, second 
and thirds integrals respectively, we get 



Mf) =f{xi)A Xl 



Xj(x I )f(x I )(x I+ i - xi) 
f(xi), 



f{xi) 



(X/+1 - Xi) 



Xj (Xi)f(xi)(xj - X/_i) (xj - X/_!) 



0{Ax 2 ) 



=f(xi)Axi - ^ ^ (xi +1 - xi) - ^ ^ f (xi - xi-i) + 0(Ax 2 ) 



--f(xi) 



Axi - ]^Axj - ^(Axj + Axi^i] 



(Ax/- Ax/_i) + e>(Ax 2 ). 



-0(Ax 2 ) 



Now we consider i = 1 

Mf) 



f(x)dx 



1/2 



Xf(x)f(x)dx. 



We apply midpoint and left rectangle rules in first and second integrals respectively 
to obtain 



Mf) =f(xi)A Xl 



Ai" (xi)/(xi)(x 2 - xi) - ^y^-(x 2 - xi) 



0(Ax 2 ) 



--f(xi) 
f(xi) 



Axi - ^(Axi + Ax 2 ) 
(3Axi - Ax 2 ). 
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Hence the proof of Lemma 13.11 is completed. □ 

Now the four main sections are described to inquire the consistency of the fixed 
pivot technique. First, we evaluate the discretization error in the birth and death 
terms in Sections 13.11 and 13.21 respectively. Then one can summarize all terms in 
section 13.31 Finally, the five different types of meshes are considered to evaluate the 
local discretization error in Section 13.41 



3.1. Discretization error in the birth term. The integrated birth term of 
aggregation PBE (|1.2I) over ith cell can be written as follows 

I rxi+i/2 r x 

Bi = -f I K(x - y,y)n(t,x - y)n{t,y)dydx. 

Let us denote 



1 [ x 

(3.2) f{t,x) = -J K(x-y,y)n(t,x-y)n(t,y)dy. 

3.1.1. Birth term on internal cells. Considering i = 2, 1 and using 
Lemma |3.1[ we can rewrite B; as 



Bi=- / / K(x - y,y)n(t,x - y)n(t,y)dydx 



\ ( x ) / K ( x ~ y,y)n(t,x - y)n(t,y)dydx +Zi{f)- 



1 

+ 2 

By changing the order of integration of the first two terms on the right hand side, we 
obtain 

i—X 



1 r*i+i/2 r*i+i 

(3.3) B,=-^ I / XJ(x)K(x-y,y)n{t,x-y)n(t,y)dxdy 

*i (x)K(x - y, y)n(t, x - y)n(t, y)dxdy 
K {x)K{x - y, y)n(t, x - y)n(t, y)dxdy 
2_ j / \ (x)K(x-y,y)n(t,x-y)n(t,y)dxdy 



j = l Jx j-l/2 
1 

2 

1 

2 

1 r x i+i/2 r x i 



x i+ i 



Xi-1/2 J Xi 
X i + i fXi+1 

+ / 



9 

j=l JX 3-l/2 



\ K (x)K(x - y, y)n(t, x - y)n(t, y)dxdy 
K (x)K(x - y, y)n(t, x - y)n(t, y)dxdy + 3i(f). 



1 

+ 2 
1 

+ 2 

Let us denote the integral terms on the right hand side in (I3.3[) by I±, . . . , I§ respec- 
tively and calculate them separately. 

The integrals Iq and I3 . First we consider the integral Iq as follows 

1 f x * 

h = 2 J 9i(t,y)dy, 

where 
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(3.4) 9i{t,y)-= \^{x)K(x-y,y)n(t,x-y)n(t,y)dx. 



Using Taylor series expansion of g(t, y) with respect to y about Xi-i, we get 

(3.5) h = gi{t,Xi-i) hg^Xi-x) — h . . . , 

4 lb 



where u. \ ®9i \ ®9i /, si 



Now we calculate t^(i,Xi_i) using the Leibniz rule to differentiate under an inte- 
gral. Then, we apply the left rectangle rule to the integrals in g[{t, Xi-i) and gi(t, x^_i) 
defined in (|3.4I) . Further, we insert A~(x^_i) = and if(0,Xj_i) = to obtain 

g i (t,x i _ 1 )=0 + O(Ax 2 ) and g'^t, x^) = + 0(Ax 2 ). 

Thus, by substituting the value of g'^t^Xi-i) and gi(t,Xi~i) in (13.5[) . we have 

(3.6) I 6 = + O(Ax 3 ). 
Similarly, we can evaluate the integral I3 as 

(3.7) I 3 =0 + O(Ax 3 ). 

The integrals I\ and I4. Let us next consider the first integral I\ from equation 
p.3j) and apply the midpoint rule to the outer integral. Furthermore we can use 
the relationship Nj(t) = n(t,Xj)Axj + 0(Ax 3 ) for the midpoint rule and substitute 
x — Xj = x' to obtain 

1 Izi nX i+1 -X j 

(3.8) I 1 = ~22N j (t) \f(x' +x j )K(x',x j )n(t,x')dx' + 0{Ax 3 ). 

j-l Jxt-Xj 

We define kj and 74 j to be those indices such that the following hold 

(3.9) Xi - Xj e A k j and := i?[(x, - Xj) - x ti J 

where TT / \ f 1 if x > 0, 

^) : =( -1 if x < . 

We will use the convention of Riemann integration that J f(x)dx = — J^" f(x)dx 
and the equation (|3.8[) can be rewritten as 



h=^Nj(t) J 



A+(a;' + Xj)K(x' , Xj)n(t, x')dx' 



i-l 'i+l,i + l(7i+i,i-l) .a: k+1/a 

X^t) 2 / a+^' + ^C^MM'M*' 



2 



^2Nj(t) / A+(x' + a; i )if(a;',a; J -)n(t,a; / )rfa;' + 0(Aa; 3 ). 



2 

J- 1 ! i + l,j+2"'i + l,3 
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Let p be the total number of terms in the following sum 

/ \f(x' + Xj)K(x',Xj)n(t,x')dx' 



fc='i,3 + !(7i,j+l) 



In particular, let p := #{n : i-y + + 1) < ft < tt-i,j + |(7i+i,j ~~ 1)} an d set 

fci := + + 1), k 2 := k x + 1, . . . , fcp_i := k x + (p - 2). 

Next, we shall show that p is finite and can be estimated by a constant which 
is independent of the grid size. By using the definition of the indices and 7,-,j in 
(|3.9j) . we can estimate 



(p - 2)Ax mm < Ax k2 + Ax ka + . . . + Ax kp _ 1 < -(Axi + Ax i+1 ) < Ax 



which implies using the assumption of quasi uniformity (|2.ip that 

(p ~ 2) < AX <C ^p<C + 2. 
^^min 

This means the above sum has uniformly bounded finite number of terms. So we 
can apply the midpoint rule to the integral in second term on the right hand side and 
use N k (t) = n(t, x k )Ax k + ©(Ax 3 ) to get 

(3.10) h=lj2Nj{t) f h ' 3 + ^ z,J X+(x' + x^Kix^x^nit^^dx' 

j = l Jxi-Xj 

1 1-1 

+ 2^2 N ^) E Xi(x k +x 3 )K(x klXj )N k (t) 

j = l Xi<Xj+X k <X i+1 

+-^JVj(t) / Xf(x' + x J )K(x',x ] )n(t,x')dx' 

^ i-1 'i+U + ^Ti + i-i^ 1 ) A 3 9 

-sE^*) E A Xl Zx t+1 0? {K ^ X ^ Xk)} 

J = 1 fe='i, 3 + |(l+7i,i) 

+C(Ax 3 ). 

Further we proceed as in /i to evaluate the fourth integral J4 on the right hand 
side from equation (|3.3[) as 

(3.11) h=y^2Nj(t) \r( x ' +x J )K(x',x j )n(t,x')dx' 

j=l Jxi^x-Xj 

1 4-2 

+ 2 E ^'W E X i ( Xk + x )K{x k , Xj )N k {t) 

j = l Xi-i<Xj+X k <Xi 



Vj2 N ^ I * X * (x' + Xj)K(x' 7 x J )n(t,x')dx / 

L Jx, , 1 



J =1 ^ij+i-^ 
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1 i-2 i i,J + i(7i,J- 1 ) 



Ax fc 3 a 



7 {if(a;fc,a;j)n(i,a;fe)} 



12^~ JV ~' ^ Axi + Axi-idx 

3 = 1 fe=i l -l,j+^(l+7.-l.j) 

+C(Aa; 3 ). 

T/ie integrals I2 and I5. To compute the integral I2 in equation (|3.3p . we apply 
the right rectangle rule to the outer integral and use Ni(t) — n(t 7 Xi)Axi + 0(Ax 3 ) 
then split it into two integrals as follows 



h =- 



I r x i+i 



Xf (x)K(x — Xi, Xi)n(t, x — Xi)Ni{t)dx 



Xi + l 



Xf{x)K(x - Xi,Xi)n(t, x - Xi)Ni(t)dx + 0{Ax 3 ). 



We use the left rectangle rule in the second integral on the right hand side and 
then put K(0,Xi) = which gives us a third order term. Therefore, we obtain by 
substituting x — Xi = x' in the first integral 



h = ijv,(t) 



X i + 1 -Xi 



X+{x' + Xi)K(x', Xi)n(t, x')dx' + 0(Ax 3 ), 



By using the indices £jj and 7^ defined in equation (|3.9p . again we proceed in the 
same way as in I\ to obtain 

(3.12) h=\Ni{i) K(xk+Xi)K(x k ,Xi)N k (t) 



Xi+X k <X i + 1 
Xi-^i —Xi 



Xf(x' + Xi)K{x', Xi)n{t, x')dx + 0{Ax 3 ). 



Before combining all these integrals together, we also need to discretize the integral 
is defined in equation (I3.3j) . Finally, analogous to the integral I2, we can evaluate is 
as 



(3.13) 



is =5^-1 (t) X l {x k +x i ^ l )K{x k ,x % - 1 )N k {t) 



Xi-l+Xk<Xi 
»Xi—Xi-i 



+C(Aa; 3 ). 



A.^ (a; + Xi-i)K(x' , Xi^i)n{t 1 x')dx' 



x, . 1 



Collecting together ([32]), dSTTJ), ffBTTU]) . ([3~TT]) . (|3~T2l) . (f3~T3l) and substituting all these 
terms into equation (|3.3[) . we obtain 



(3.14) 



B, 



Y%L ■, . , . (l-^ k yt(x k +x 3 )K(x k , Xj )N 3 (t)N k (t) 
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+ \ E -=i Nj(t) SvX™ W + Xj )K(x>, Xj )n(t, x')dx' 



E 



+2 E 5=i W) J^T-C 2 ^ 1 '* x ^ x ' + X ^ K ( X '> x '^ dx ' 

+ \ E}=i (*) C~ X \ K( x ' + Xj)K(x', Xj )n(t, x')dx> 



12 Z^,=i jv jW Z^ fc=ii 3 + i (1 



t,i + *T<, 



i*+i,i + 5(7»+i,3'— 1) Aa; fe 



Aii+Ai i+ i 3a: 

+0(Ax 3 )+Uf)- 



The first two terms on the right hand side are exactly the fixed pivot discretization 
Bi from (|2.9[) . Let us denote the sum from third to sixth terms and the difference 
of remaining two terms on the right hand side by E and E' , respectively. Doing 
so and substituting Zi(f) from Lemma 13.11 equation Q3.14p can be written as for 
i = 2,. 



(3.15) B l =B l + E + E 

f(Xi) 



A.I 



V 2 



12 



(Ax l+ i - Axc*_iK Axi 



Axi-i + Ax i+1 



0(Ax 3 



3.1.2. Birth term on boundary cells. Now we evaluate the order of the 
error in birth term on boundary cells. We will treat the boundary cells separately. 
Therefore, first we consider the birth term for % = 1. Using the Lemma 13.11 and 
changing the order of integration we get 



flx=i 

2 



\^{x)K(x - y, y)n(t, x - y)n(t, y)dxdy 



Xi/2 J X t 

I r x i+i/2 r^2 



K [x)K{x - y, y)n(t, x - y)n(t, y)dxdy 



Xl 

xi rX2 



1+1/2 J V 



\f(x)K(x - y, y)n(t, x - y)n(t, y)dxdy + Zi (/)• 



We apply the right rectangle rule to the outer and inner integrals respectively of the 
first term as well as to the outer integral of the third term. Both of these terms are 
of second order separately. Also we use the left rectangle rule to the outer integral of 
the second term to get 

1 f X2 

Bx = -Axi / \t(x)K(x-x 1 ,x 1 )n(t,x-x 1 )n(t 7 xi)dx + 0{Ax 2 )+Zi{f)- 

4 J xi 



Again by applying the right rectangle rule and putting N\{£) = n(Jt,x\)Axx 
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0(Ax 3 ), we obtain 

Bi = ijVi(t) X ^ Xk +xi)K(x k , Xl )N k (t) +0{Ax 2 ) +Zi(f). 

Xl+X k <X2 

In terms of fixed pivot discretization we have 

(3.16) B 1= Bi + ^|^[3A^i - Ax 2 ] + 0(Ax 2 ), 

where / is defined in Eq. (|3.2j) . Then the application of the right rectangle rule to 
the integral in (|3.2j) gives f(xi) — O(Ax). Substituting this value of f{x\) in (|3.16|) . 
we obtain 

(3.17) B l= B l +0{Ax 2 ). 

Finally, in case of last boundary cell i.e. i = I, we use the Lemma l3.il and change 
the order of integration as in B\. Further we split each integral into two parts to 
estimate 

i-i 



Bi =-^2 / Xj(x)K(x-y,y)n(t,x-y)n(t,y)dxdy 

Xj(x)K(x - y, y)n(t, x - y)n(t, y)dxdy 



k=l Jx k-l/2 •> XI 
Y f x I f x I+l 



'Xl-X/2 Jx I 

+- / Xf(x)K(x-y,y)n(t,x-y)n(t,y)dxdy 

* J xi Jy 

^ pxi+i r x i+i 



XJ (x)K(x - y, y)n(t, x - y)n(t, y)dxdy 

x I+l/2 J V 

Y [~~ 2 r x k+i/2 r x i 

+ o^2 / Xj{x)K(x-y,y)n(t,x-y)n(t,y)dxdy 

k=l J x k-l/2 Jxi-1 

^ rxi-i pxi 

+- / Xj(x)K(x-y,y)n(t,x-y)n(t,y)dxdy 

Jx I _ 3/2 Jxi-l 

1 r x i~i/2 r x i 

+ - / XJ (x)K(x - y,y)n(t,x - y)n(t,y)dxdy 

z Jxi-i Jy 

2 r x i r x i 

+- / XJ(x)K( X -y,y)n(t, X -y)n(t,y)dxdy + 3i(f). 

J x i-i/2 J y 

By applying the midpoint, right, left and right rectangle rules to the outer integral 
of the first and fifth, second and sixth, third and seventh as well as fourth and eighth 
terms respectively. Then after substituting Nk(t) — n(t, Xk)Axk+0(Ax 3 ), we obtain 



1 i~ \ r Xi+i 
1 k=i Jx < 



(x)K(x — Xk, Xk)n{t, x — Xk)dx 



+—Nj(t) / Xj(x)K(x — x/,xi)n(t,x — xi)dx 

* J XI 

1-2 ,X f 

[x)K(x - Xk,Xk)n(t, x - Xk)dx 



■\j2N k (t) r x 7 

z k=i Jxi - 1 
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+-JV>_i(*) / Xj (x)K(x-x I - 1 ,xi- 1 )n(t,x-x I ^ 1 )dx + 0(Ax 2 )+Z I (f). 

The first, second, third and fourth term on the right hand side can be solved 
similar to I\, -Z2, I a and 1$ respectively. Thus, by substituting the value of Zi(f) from 
Lemma |3. 11 we have 



(3.18) 



A = A 



/(*/) 



[Ax j - Axr_i] + 0{Ax 2 ). 



3.2. Discretization error in the death term. From equation (J23J), the inte- 
grated death term can be written as follows 



A = 



E 



X3+1/2 



K(x, y)n(t, y)n(t, x)dydx. 



The application of the midpoint rule to the outer and inner integrals gives us 



(3.19) A = Niitf^KixuxrfNjty) + 0(Ax 3 ) = D z + 0{Ax 3 ). 

j=i 

3.3. Summary of all terms. Finally from the equations (|3.15[) . (I3.17p . (|3 . 1 8|) 
and (|3.19p . we can summarize the spatial truncation error <7j(t) = A — A — (A — A) 
as follows 



(3.20) 
(3.21) 



(3.22) 



01 (t) = G(Ax 2 ), 



o-i(i) =E + E 
f'{x. 



12 



Ax t - 



(Ax i+ i - Axi^i)< Ax 



Axi-! + Ax i+1 



Axt-! + Ax i+1 



+0(Ax 3 ), i = 2,...,I-l, 



ai(t) 



[Axi - Aa)/_i] +C(Aa 



where -E and A are defined in p,14j) . 

3.4. Meshes. Now let us consider the following five different types of meshes to 
calculate the order of local discretization error. Most of the grids are the same which 
are used for solving breakage PBE in J. Kumar and Warnecke [24] . 

3.4.1. Uniform mesh. Let us begin with the case of uniform mesh i.e. Axi — 
Ax and Xi — (i — l/2)Ax for any i = 1, To estimate <7j(t), first we have to 
evaluate the order of E and E'. In case of such uniform grids, we have 



-i+l/2) x i+i ~ x j 



c i-j+3/2) 



and 



Xi—l Xj — Xi_j_l/2 • 



By using the definition of indices in ()3.9j) . we have 



c i-j + l/2 



e A, 
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This implies that Xi - x 3 = Xi- j+1/2 = x h j + 1 7 . . . 

Similarly, we can easily obtain 

Xi+i - Xj = x l _ J+3/2 = x u+i ] + p 

aild " X 3 = = X U^+^ t - ld - 

Therefore, from Q3.14p we have E = 0. For uniform grids, we can observe from above 
that and Xi+i — Xj are the right end boundaries of adjacent cells, 

i.e. h-ij = (i — j — l)th, kj = (i — j)th and k+ij = (i — j + l)th cells, respectively. 
Here Ji-i,j = = 7»+ij = 1- Thus, by substituting the values of all these indices 
in E' defined in (|3.14[) and using the Taylor series expansion we obtain E' = ©(Ax 3 ). 
Finally, from the equations ([3T20)l - (f3T22]l we estimate 



0(Ax 2 ) i=l,I, 
0(Ax 3 ) i = 2,...,I 



The order of consistency is given by 

lk(<)ll = M*)l + E + = °( A ^ 2 )- 

i=2 

Thus the method is second order consistent on a uniform mesh. 

3.4.2. Non-uniform smooth mesh (2^+1/2 = rx i _ 1 / 2 , r > 1, i = 1, . . . , I). In 

general, non-uniform smooth grids can be obtained by applying some smooth trans- 
formation to uniform grids. Let us consider a variable £ with uniform grids and a 
smooth transformation x = g(£) such that ^±1/2 = g{d±i/2) for & n y * = 1, . . . , / to 
get non-uniform smooth mesh. The scheme again gives us second order accuracy. Let 
h be the uniform mesh width in the variable £. In case of smooth grids, Taylor series 
expansions in smooth transformations give 

A Xi = x l+1/2 - Xi _ 1/2 = .9(6 + h/2) - gfa - h/2) = hg'fa) + ^g"&) + 0(h A ). 

Hence, by evaluating Axi-\ and Axi+\ in a similar way, we obtain for any i = 2, . . . , I 
Axi - Ax^i = 0(h 2 ), 2Axi - (Axi-i + Ax i+1 ) and Ax 2 - Ax^Ax.^i = 0(h 3 ) 
The above identities help us to simplify the equations (|3.20[) - (I3.22|) and give 

(3.23) a 1 (t) = 0(h 2 ), a I {t)=0(h 2 ) and a^t) = E + E' + 0{h 3 ). 

To find the order of <7j, we need to evaluate E and E' separately. Let us first consider 
E from (I3.14[) and set f(t,x,y) := K(x,y)n(t,x). Then the applications of the left 
rectangle rule to the integrals in first and third terms as well as the right rectangle 
rule to the integrals in second and fourth terms estimate 

1 x — (^/- +--V- ■ Xi ~\- Xj} 

(3.24) E=-J2 Nj(t) " L — /(*, ^ - x h x 3 ) 

._y Xi Xi-\-\ 

\ (xi ii-, — Xi+1+Xj) 2 
~a l^ N iW f{t,x i+1 -Xj,Xj) 



16 



ANKIK KUMAR GIRI AND ERIKA HAUSENBLAS 



1 i-2 

zEw*) 



i-l 



3=1 
i-l 



-«EW) 



3=1 



■E^w 



3=1 
i-2 



-sE^w 



3=1 
i-l 



Xi 



Xi Xi — 1 



^z+l 



#2 



/x' X{ Xj , ) 



6 



3=1 



^i Xi — \ 



f X '(t,Xi - Xj,Xj) + 0{Ax 3 ). 



Let us denote the combination of first four and last four sums on the right-hand side 
in p. 241) by E\ and E% respectively. We shall solve them separately. First we take 
Ei and approximate / at (t,Xi — Xj,Xj) by / expanded around (t, x^i — Xj,Xj) in 
first and fourth summations to get 



£i = 2 ][>•(*) 



3=1 



Xi ^i+1 



X{f(t,x l+1 - Xj,Xj) + (Xi - X i+ l)f x t(t,X i+1 - Xj,Xj)} 



I - ■ • K +M +| 7 »+i, 3 ' " ^+i + x if 



■E^w- 



3=1 
i-2 



*£i ^i+1 



f(t : Xi-\-i Xj : Xj) 



+zE^w 



3=1 
i-l 



-zE^-w 



3=1 



*£i X^— \ 



f(t,Xi-\ Xj,Xj) 



:{f(t,x t -i - Xj,Xj) + (xi - Xi-i)f x >(t,Xi-i - Xj,Xj)} + 0(Ax 3 ). 



Again approximating f x i at point (t, Xi+x — Xj,Xj) by f x > evaluated at point 
(t, Xi-i — Xj, Xj) in the second part of the first summation. Then we replace j by j + 1 
and j — 1 respectively in second and third term summations. Further we can use the 
relationship Nj(t) — n(t, Xj)Axj + 0(Ax 3 ) for the midpoint rule to obtain 

E x =~n(t, xi)Axi f(t, Xi+i -xx,xi) 

Q Xi-\-\ Xi 



i 



zE n ^ 

x j+ i)Ax, 



3+1" 



3=1 
i-l 



K +1 , J + 1 + l 7l+1 , J + 1 - Xi+l + Xj+l)' 
X i+ 1 - Xi 



■f(t,X i+ i - Xj+l,Xj+l) 



E(Xi^ .j r L^ i ■ Xi -\- Xj^j 
n(t,Xj)Axj — - 2 ^ ' ,J f(t,x i+1 -Xj,Xj) 



3=1 



Xi-\-l Xi 
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1 1-4 ( X h-i ;_i + ± 7i _i ;_i ~~ X %-1 + X 3-l) 2 

+ -2_^n(t,Xj-i)Axj-i — - — — — — '— f(t,Xi-i - Xj-i,xj-i) 

J — z 

1 i X l- - + -7- ■ — X i X j)^ 

-j2^n(t,Xj)Axj — 1,3 2 ^- f{t,Xi-i - xj,Xj) 

1 f 3// 1 ■ X " "I - X\ ") 

-±n(t, a^As / ' M+a7M -/(t, - an, si) + G(Ax 3 ). 



Approximating the functions x' i-» f(t, )n(i, a;') at point Xj by /(t, arj±i — 

x',x')n(t, x') evaluated at point x' = Xj±\ of the third and fifth terms on the right 
hand side respectively. Then we approximate / at (t, Xi-i — Xi,xi) by / expanded 
around (t, Xi+\ — x\,x{) in the last term on the right-hand side and rearrange the 
whole expression for E\ as follows 



Ei = ^-n(t, xi)Ax ^ l ' t '' 1+i1 ' +1 ' — f(t,Xi+i - xi,xi) 



4 Xi-^-i Xi 

1 1-1 

+ 7 ^2 n (^ x i+^( Ax 3+^ - Ax j) 



3=1 



X ^ 2 ^_ f(t,X i+1 -X j+1 ,X j+1 ) 

Xi-\-l Xi 



J = l 

1 i_1 d 

{x Xj . . 2 

(X/. . 1 i_, . X{ -r Xw) 



i-1 



1 

K-U-i + ^i-Li-i -Xi-l+Xj-l) 



j (tj X^_i Xj — i , Xj—X ) 
x^ — x^—l 



X 

1 i_1 a 

i=2 



(^/ . ■4-l'v. ■ X t -\- X*,) 

<r. -j, . H,j -r 2 it, 3 J ' 
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-jn(t, x 1 )Ax 1 ^ h - 1+h ^ Xl+Xl) f(t, x i+1 -x u x{) + 0(Ax 3 ). 

4t Xi Xi — \ 

Let us denote the each term on the right hand side by En,...,Eig respectively. 
Therefore, we have 

(3.25) E 1 = {E n - E ia ) + (E 12 - E 15 ) + (E 13 - E 16 ) + (E u - E 17 ) + 0(Ax 3 ). 

In order to solve the error E\ , we need to evaluate differences of the above mentioned 
terms separately. First of all, we are interested in £12 — -E15. In case of smooth grids, 
we have as before Axj — Axj-i — 0(Ax 2 ) which clearly shows that £15 = 0(Ax 2 ). 
Then Taylor's series expansion gives us a second order approximation of £15 as 



1 4-1 

E15 = j}n(xj + i)f(xi + i - x j+ i,x j+ i)(Axj - Axj-i) 



J=2 

x V l ^zl±l2izM^ Li L^L + (Ax 3 ). 

X% X% — J 

Without loss of generality the summation in E 12 can be started from j — 2 since the 
term we are dropping from the summation is third order accurate. Then we can write 



(3.26) E u - E 15 = n(xj +1 )f(x i+1 -x j+1 ,x J+1 ) 



3=2 



x 



f (Ax i+1 - Axj) 2 
1 Ax l+ Ax l+1 + 

+C(Aa; 3 ). 

By applying Taylor series expansion in the smooth transformation, we obtain 

(3.27) (Axi + Aasi_i)(Aa: i+ i - Aay) = 2h 3 g'(^)g"^ j ) + 0(h 4 ), 
and 

(3.28) (Aa^ + Ax l+1 )(A^- - Aa^-i) = 2/i 3 </ (&)<?"(&) + <^ 4 )- 

Here we consider a particular type of non-uniform smooth grids i.e. 2^+1/2 = nEj_i/2, r > 
l,i = 1, ...J. Such grids are called geometric grids. These grids can be ob- 
tained by applying a smooth transformation as exponential function on uniform 
grids. Here we have x i+1/2 = exp(& +1/2 ) = cxp(/i + Ci-1/2) = exp(/i) exp(^_ 1/2 ) = 
exp(/i)x i _ 1/ 2 =: rXi_i/ 2 , r > 1. 

Let us consider £ n , £ 12 , £217 £22, £31 an d £32 are corresponding points on uniform 

mesh for x h+i> ^+1 - ^+1, ^, 3 + i 7liJ - " ^, !B l 4 - 1 j- 1 +i7*-ij- 1 and 
— Xj-i, respectively. These are defined as follows. 

61 =ln(^ i+1 , J+1 +i 7<+1<3 . +1 ).---.C32 = lnfo-i - arj-_i). 

By the definition of the indices in (|3.9[) . we know 

Xj+i - x j+1 G A Ji+1J+1 , ^ - ccj G A(! 4ij and Xi-i - xj-i G Ai i _ li;f _ 1 . 
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For geometric grids, we have 

x l+1 - x j+1 = r(xi - Xj) = r 2 (xi-i - Xj-i). 

Therefore, we have 

h+l,j+l = h,j + 1 = h-lj-l + 2. 

Further, in case of geometric grids, we have 

7»+i,i+i = 7»,j = Ti-ij'-i- 

Let us consider 



hi = - £12 =H x i i+1 . 3+1 +i 7l+1 , J+1 ) - Hxi+i - Xj+l) 

Xi+i - Xj+i J \ Xi — x 

.-i.j-i+^7i-i. 
Xi—1 Xj—i 



In! -— — — ^— — - I 



Similarly, we can estimate 

£ia - 62 = Inf 3 ^ 1 " 3 ^ 1 ) = ln(r 2 ) = 2h. 

\Xi-l - Xj-! J 

By the application of smooth transformation we can write 

(3-29) [x h+l j+1 + H+1]+1 ~{x l+1 - x, ]+l )f 

= Wxi) .9(£i 2 )] 2 = h\{g'{i 12 )Y + 0(h 3 ), 

and 

( 3 ' 3 °) ^I-l.-l + iT.-i,,-!"^- 1 ~ Zj-l)] 2 

= [.9(6i) - 5fe 2 )] 2 - h\{g>(&)} 2 + 0(h 3 ). 

Collecting (EOT)) . (jX28l) . (pH9)) and (f3T30l) together, we obtain 

(Ax i +Ax i - 1 )(Ax j+1 - Axj)(x h+lt . +1 + i li+1:j+1 - x i+ i + x J+1 ) 2 

-(Ax, + Ax l+ i)(A Xj - Aa; j _i)(a; i ._ 1(J _ 1+ i 7 ._ li - + ^--i) 2 
=-8/ l 4 / llff 'te)5"fe)3'(a 2 ).9"(a2) + O^ 6 ) = 0(^ 6 ). 

Therefore, by substituting the preceding identity in (|3.26p we estimate 

(3.31) E 12 -E 15 =0( Ax 3 ). 

Next, we need to compute E13 — £7x6 defined in (|3.25|) . In case of smooth grids, 
by using smooth transformations and Taylor's series expansions, it is easy to show 
that 

k J+ i 7M - ( x i - x i)\ 2 - fai-i.i-i+^i-i.i-i - (^-1 - x j-i)f = 0(h 3 ). 
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which implies that i5 16 = 0(Ax 2 ). So we can take a second order approximation of 
£16 by using Taylor's series expansion as 



^-1 



fjXj+l - x j+ i,x j+1 ) 
Axi + Axi-i 



X <^ (^/ iii + i 7iJ ~ Xt + Xj) - {xi^ ^, 



J'=2 



+ 27i-i,i-i 



Xj,— 1 "T X j—\ y 



C(Ax 3 ) 



Similarly as before, we can have 

Therefore, again without loss of generality the summation in E13 can be started from 
j = 2 since the term we are dropping from the summation is third order accurate. 
Let us estimate 



(3.32) 



E13—E16 

i— 1 



1 



j=2 



Azi + Aa; i+ i 

+ Xj) 2 - (X l ^ 1 + 1 - !Ei_l + Xj-i) 2 } 



A^i + A^-i 
+C(Aa; 3 ). 

By using smooth transformation and Taylor's series expansions, we can easily obtain 

(A.i 1 +Ai ! _i){( 3 ; 1i+1j+i + i 7i+1j+i - x l+1 + x 0+ if - {x lij+ i^ 3 -Xi+ Xj) 2 } 

-{Axi + Aa5i + i){(a; j4i;(+ i 7ij;j - x t + Xj f - (a; Ji _ liJ _ 1+ i 74 _ 1|;f _ I - Zi-i + ^-i) 2 } 
=2h 2 h\g>{i 22 )g"{i 22 ) - 2h 2 h 2 g'^ 22 )g"^ 22 ) + 0(h 5 ) = 0(h 5 ). 

Putting the above identity in equation Q3.32p . we have 

(3.33) E 13 - E 16 = 0(Ax 3 ). 

Now we shall approach to evaluate E14 — En. Since the terms appearing in the 
summation En are third order accurate. Therefore, the summation in £14 can be 
started from j = 2. By using Taylor's series expansion in En we consider 

1 i_1 d 

(3.34) E u - E 17 = -S2Ax 3 —{n(x J+1 )f(x l+1 - x ]+ll x J+1 )} 



3=2 



dx' 

Axj + Axj + i Axj + Axj-i 



Axi + Axi+i Axi + Axi_i 
+C(Aa; 3 ). 

Analogous to equations (|3.29p and (|3.30p . we estimate 



(3.35) 



x. + xjf = 0{h 2 ). 
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By applying the smooth transformation and Taylor's series expansion, we have 

(3.36) (A Xl + Axi-x^Axj + Ax j+1 ) = 4ftV(6)fl'(^) + 0{h 3 ), 
and 

(3.37) (Axi + Ax i+1 )(A Xj + Ax^x) = 4/iV(&)</(^) + 0(h 3 ). 
Substituting the identities (|3.35[) , (|3.36p and (|3.37[) in equation (|3.34[) , we obtain 

(3.38) E 14: -E 17 = 0(Ax 3 ). 

Finally, we have to estimate En — E\s in equation (|3. 251) . Since the geometric grids 
considered here are monotonically increasing. So, the points Xi — x\ and aJi+i — x\ 
lie in ith and (i + l)th cells, respectively. By using the definition of the indices lij 
and in (|3.9p . we have i^i = i, U+i.i = i + 1 and 7^1 = 7i+i,i = — 1. Then the 
expression for En — E\g can be written as 



(3.39) E n -E m = -n(t,x 1 )Ax 1 



(Xi+l/2 - + Xif 



Axi + Ax i+ i 

A 2 1 

f(t,X l+ i - Xi,Xi). 



{ x i-l/2 ~ x i + x l) 2 



Axi + Axi-1 

Let us use the following formula and the identities mentioned before (|3.23[) to get 
(Axt + A Xl -i)(x i+1/2 - x i+ i + xi) 2 - (Axi + Ax l+1 )( Xl _ 1/2 - x { + xi) 2 

=i[Axj(Ax i+ i + Axi)(Ax i+1 - Ax % ) - Axi 2 (Aa; i+ i - Aa;»_i) 
+ (2Axi - Ax l+1 )(Ax 2 - Axi^iAxi+i)] 

Inserting the above identity in equation (|3.39|1 . we estimate 
(3.40) £n - E 1S = 0(Ax 3 ). 

Finally, we substitute (j3~3T|) . $£55$ , P~351) and (pT4H)) into (j3~25)l to obtain 

E x = 0(Ax 3 ). 
Similar to E\, we can evaluate E2 defined in (|3.24[) as 

E 2 =0(Ax 3 ). 
Therefore, substituting E\ and E2 in (|3.24[) . we obtain 

E = 0(Ax 3 ). 

Now let us consider E' from (|3.14|) and set f x {t,x,y) :— -^{K(x,y)n(t,x)}. Then 
we replace j by j — 1 in the second term on the right hand side and apply Taylor 
series expansion about xj = Xj-i in n(t, Xj)f x i{t, Xk, Xj) of the first term. Further we 



22 



ANKIK KUMAR GIRI AND ERIKA HAUSENBLAS 



rearrange the expression for E 1 to estimate 

^ i-i h+i.j+hhi+ij- 1 ) 
E'=— ^n&x^Axj^-Axj) ^ 



Ax k d 



3=2 
i-1 

12 



fc=«i,i+i(i+7i,i) 



Ax* + Ax i+ i 



fx'{ti X kl Xj — i ) 



ii+l,j + ^(7i+l,j-l) 

L2 ^ J ^ Ax, + Ai !+ i 

J=2 fc=Jij+i(i-Hr«,i) 

i-1 'i,j-i + 5(7»,3-i-l) 

-Y,n(t,Xj-i)&Xj-i E 



f x ' (t, Xfc , Xj — i ) 



Ax fe - 



12 



fc=Zi-l,j-l + i(l+7i-l,j_l) 



Ax, + Axi_i 



/a;' (t> %k: l) 



J=2 

+C(Ax 3 ). 

For smooth grids, we have Axj_i — Ax.,- = 0(Ax 2 ). Therefore, the first term on the 
right hand side of the above equation is of third order. By using the Taylor series 
expansion about Xk = Xj in f x > (i, x kl Xj-i) in the remaining terms on the right hand 
side, the above equation can be further rewritten as 

'i, 3-1+3(7*, j-i-i) 



(3.41) 



1 

E' = — ^n(i,^_i)A^_i 



i=2 



E 



Ax,,-' 



h+i,j + i(li+i,i-l) 

E 

fc=i,, 3 + i(l+ 7 ,, J ) 



fe=/»-l,3-l + 2 (l+7»-l,3-l) 

Ax k 3 



Axi + Ax i+1 



Axi + Ax^i 
/^(t,x J ,x,_ 1 ) + 0(Aa; 3 ). 



For non-uniform smooth grids of the type x i+1 / 2 = r%i-i/2> r > 1, i = 1, ... ,1, we 
have 

x i+ i - xj = r(xi - xj-i) and x i+ i - Xj e A; i+1 ^, ij-ij-igA[ (iH . 

Therefore, we obtain Zj+ij = Zi,j-i + 1- Similarly, we can have Zjj = Zi-i.j-i + 1. In 
case of such smooth grids, it can easily be seen that 7i-ij-i = 7,,j and Ji,j-i = 7»+i,j- 
Let fci := + |(1 + 7i_i,j_i), fc 2 := + |(1 + 7»,j) = fci + 1, . . . , fc p := 

+ 2"(7ij-l — 1) = + 1; : = + 2^ i + 1 'j 

can write 



1) = k p + 1. Then we 



(3.42) 



«i,i-i + 5(7M-i-l) 

E 



Ax fc " 



i i+i,j + 3(7i+i,j- 1 ) 

E 



Ax; + Axi_i ^—^ Axi + Ax i+ i 

fc=I ( -ij-i + i(l+7«-ij-i) *='ij + i(l+7*,i) 



Ax fcl : 



Ax fc2 



A.t 4 + Axj_i Ax* + Aa; i+ i 
Ax fe 3 



+ ... 



+ 



Ax fcp+1 : 



= E 



Ax*. 



Axi + Axi-i Axi + Ax i+ i 



Axu 



Axi + Axj_i Axi + Axi+i 



where fc m +i = k m + 1. Let £ m i, £ m 2, £ m 3 are the corresponding points on the uniform 
mesh for x km+3 / 2 , x km +i/2 an d x km _ 1 / 2 , respectively. Since 



In 



x k m +3/2 
x k m + l/2 



ln(r) = h = In 



£fc m + l/2 
X fe m -l/2 



m2 sm3 • 
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Then by using smooth transformations and Taylor's series expansions, we can have 
(3.43) (Axi + Ax l+1 )Ax 3 km - {Ax, + Axi-i)Ax 3 km+1 = 0(h 5 ). 

Therefore, by using (|3.43[) and (|3.42l) in equation (|3.41l) . we estimate 

E' = 0{Ax 3 ). 

Thus, we insert the value of E and E' in equation (|3.23|) to obtain 

a, = 0(h 3 ). 

Finally, analogous to the uniform mesh, the technique is second order consistent. 

3.4.3. Locally uniform mesh. Let us explain an example of a locally uniform 
mesh. First, the computational domain is split into many finite sub-domains and each 
sub-domain is further split into an equal size mesh. This gives us a locally uniform 
mesh. It is not easy to analyze the order of consistency on locally uniform mesh. 
So we calculate it later numerically. In this case, the scheme gives only first order 
consistency. 

3.4.4. Oscillatory mesh. A mesh is called oscillatory mesh, if for any r ^ 1 > 

0, we have 

. f r Axi if i is odd, 

l\Xi + i := < i . . . 

[ -Axi if i is even. 

Since there is no cancellation in the leading error terms of equations (|3.20p - (|3.22p 
as well as£ = O(Ax) and E' = 0(Ax 2 ) we have 

tr<(t) = 0(Ax), i = l,...,I, and \\a(t)\\ = 0(1). 

Therefore the fixed pivot method is unfortunately inconsistent on oscillatory meshes. 

3.4.5. Non-uniform random mesh. Finally the scheme is examined on non- 
uniform random grids. Similar to the case of oscillatory mesh, we have ||cr(f)|| = 0(1). 
Thus the method is again inconsistent on non-uniform random meshes. 

4. Lipschitz conditions on B(N(i)) and D(N(t)). Let us consider the birth 
term for < t < T and for all N, N e M. 1 

i 

||B(N(t)) - B(N(t))|| = ]T mm)) A(N(t))|. 
<=1 

From (|2.10l) . there exists a L > such that K(x,y) < L for all x, y €]0,a;max]- 
Then by using this upper bound L and < Xf(x) < 1 from the definition in (|2.8|) . 
we obtain from (|2.9p 

||B(N(t))-B(N(t))|| 4 L EE E mQNkW-NjWMW 

1=1 j — 1 Xi<.Xj-\-Xk<.XiJ r x 

+ ^ L EE E W^N^-N^N^)] 

1=1 j=l X i - 1 <Xj+X k <X i 

I I 

j=l k=l 
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Now we enjoy a useful equality Nj (t)N k (t)-Nj {t)N k (t) = \ [(N } (t)+Nj (t))(N k (t)- 
N k {t)) + (Njit) - NjityiNkit) + N k (t))] to get 



(4.1) 



|B(N(t))-B(N(t))||<ii^^ 



j=l k=l 



\(N 3 (t) + N 3 (t))\\(N k (t) - N k (t))\ 
+\(N j (t)-N j (t))\\(N k (t)+N k (t))\ 



It can be easily shown that the total number of particles decreases in a coagulation 
process, i.e. Y^j=i -^Y?'(*) — := Total number of particles which are taken initially. 
Therefore, the equation (|4.1j) can be rewritten as 



(4.2) ||B(N(t))-B(N(t))||<i\#L 



- Mm + E iww - ^(*))i 

k=l j = l 

< 2N%L\\N{t) - N(t)\\. 



Finally, we consider the death term 



|D(N(t)) - D(N(t))||= E lA(N(t)) - A(N(*))| 



i=l 3=1 

t=i j=i 

Again we use the same equality as before to get 

(4.3) ||D(N(t)) - D(N(i))|| < 2iV°£||N(i) - N(t)||. 

Therefore, the application of Theorem 12.41 implies the convergence of the fixed pivot 
technique for aggregation PBE (11.21) and the convergence is of same order as the 
consistency. 

5. Numerical examples. We now justify our mathematical results on the con- 
vergence by taking a few numerical examples where we numerically evaluate the ex- 
perimental order of convergence (EOC). A detailed comparison of numerical results 
of number density and moments with analytical solutions can be found in |25l 121] . 
All the numerical experiments are performed as in |24j . 

First, we consider the following normally distributed initial condition (NIC) 



n(0, x) 



1 



: exp 



2a 2 



OV27T 

In addition, we take the following aggregation sum and product kernels 
(5.1) K(x,y) = k (x + y) and K(x,y) = k xy. 
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Since analytical solutions are not available for the above initial condition and aggre- 
gation kernels, we use the following formula in order to calculate the experimental 
order of convergence 

EOC = mf y M2 ). 

\\\N h/2 -N h/4 \\J / 

Here N/j represents the numerical solution on a uniform mesh of width h. The other 
parameters are a 2 — 0.01, fi = 1 and ko = 1. Now we will calculate the EOC on five 
different types of uniform and non-uniform meshes. 

Let us first calculate the EOC for uniform meshes. For a uniform mesh, we fix 
the minimum and maximum values of x as and 15, respectively in the numerical 
computation. The number of grid points is denoted by GP in the following tables. 
The numerical results are shown in Table 15.11 As expected from the mathematical 
analysis, the numerical results show the convergence of second order. 

Let us now consider the second case of non-uniform smooth meshes. In particular, 
we took the case of geometric grids which can be obtained by applying an exponential 
smooth transformation as x — exp(£). Here £ is the variable with uniform grids. The 
computational domain in this case is set as [le — 6, 1000] which corresponds to the £ 
domain [ln(le — 6), ln(1000)]. The numerical results have been summarized in Table 
15.21 Once again the numerical results show that the fixed pivot technique gives second 
order convergence on non-uniform smooth meshes. 



Table 5.1 Table 5.2 

Uniform grids (NIC) Non-uniform smooth grids (NIC) 



(a) K(x,y) = 


k (x + y) 


(b) K(x,y) 


= k xy 


(a) K(x,y) = 


k (x + y) 


(b) K(x, 


y) = k xy 


GP Error L 


i EOC 


Error Li 


EOC 


GP 


Error L 


i EOC 


Error Li 


EOC 


60 








60 










120 0.0598 




0.0306 




120 


0.0456 




0.0374 




240 0.0178 


1.75 


8.4E-3 


1.86 


240 


0.0118 


1.95 


9.6E-3 


1.96 


480 5.0E-3 


1.82 


2.3E-3 


1.89 


480 


3.0E-3 


1.97 


2.4E-3 


2.00 


960 1.3E-3 


1.95 


6.0E-4 


1.95 


960 


7.6E-4 


1.98 


6.0E-4 


2.00 



The third test case has been performed on a locally uniform mesh using the 
same computational domain as is in the previous case. In this case we started the 
computation on 30 geometric mesh points, and then each cell was divided into two 
equal parts in the further refined levels of computation. In this way we obtained a 
locally uniform mesh. The EOC has been summarized in Table l5~3l Table l5~3l clearly 
shows that the fixed pivot technique is only first order accurate. 

Now we consider the fourth case of an oscillatory mesh to evaluate the EOC. Let 
us take an example of oscillatory mesh, i.e. 

. f 2Axi if i is odd, 

l\Xi + i := < i . .„ . . 

I ^i\Xi it i is even. 

Here the computational domain is the same as for the first case. First, we divide 
the computational domain in 30 equidistant mesh points, and then each cell into two 
parts with 1 : 2 as per further refined levels of computation. The numerical results 
has been shown in Table 15.41 As expected, Table 15.41 exhibits that the fixed pivot 
technique is not convergent on oscillatory meshes. 

Finally we consider the fifth case of a non-uniform random mesh. The computations 
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Table 5.3 
Locally uniform grids (NIC) 



(a) K(x,y) = k (x + y) (b) K(x, y) = koxy 



GP 


Error Li 


EOC 


Error L\ 


EOC 


60 










120 


0.0416 




0.0254 




240 


0.0212 


0.97 


0.0126 


1.01 


480 


0.0105 


1.01 


6.0E-3 


1.07 


960 


5.1E-3 


1.04 


2.8E-3 


1.09 
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Table 5.4 
Oscillatory grids (NIC) 



(a) K(x,y) = k (x + y) (b) K(x,y) = k xy 



GP 


Error L\ 


EOC 


Error L\ 


EOC 


60 










120 


0.0650 




0.0347 




240 


0.0632 


0.04 


0.0309 


0.16 


480 


0.0523 


0.27 


0.0257 


0.26 


960 


0.0518 


0.01 


0.0253 


0.01 



have been performed on the same domain as for the second case. We started again 
with 30 geometric mesh points, and then each cell was divided into two parts of 
random width in the further refined levels of computation. For each value of / = 
60, 120, 240, 480, we performed five runs on different random grids and the relative L\ 
errors were measured. The mean of these errors over five runs is used to calculate the 
EOC. The numerical results have been shown in Table [5~5l Table [5751 shows clearly 
that the fixed pivot technique is not convergent. 

Next, we take an exponentially decreasing initial condition (EIC), namely n(0, x) = 
exp(— ax). The sum and product aggregation kernels in (|5.1[) are again considered. 
Since the analytical solution is known for the above initial conditions and kernels 
which can be found in |37[ [T] . Then the experimental order of convergence can be 
determined by the formula EOC = \n(Ej / E21)/ ln(2), where Ej and E21 are the L\ 
error norms. The subscripts / and 21 correspond to the degrees of freedom. We can 
calculate the error Er on a mesh with / cells. The relative error has been calculated 
by dividing the error ||N — N|| by ||N||. The parameter a = 10 is taken in the above 
initial condition. Again we computed the EOC on five different type of meshes as 
before. 

In case of uniform mesh, we set the computational domain as [0,30] to evaluate 
the EOC numerically. The numerical results are summarized in Table [5761 Again, we 
obtain the convergence of second order numerically. 

Table 5.6 
Uniform grids (EIC) 

(a) K(x,y) = k (x + y) (b) K(x,y) = k xy (a) K[ ^ y) = k(){x + y) (b) K ^ y) = k(jXy 
GP Error L x EOC Error L, EOC GP Error L\ EOC Error Zj EOC ' 

120 0229 - 0112 - 60 °- 0486 " °-° 274 

240 lolrl -0.71 9T 0.25 ^0 0.0135 1.84 7.2E-3 1.92 

480 0.0406 -0.11 0.0135 -0.51 ^4 . E 3 .94 . E-3 .92 

960 0.0402 0.01 0.0129 0.06 — ±}^Jt 



Table 5.5 
Non-uniform random grids (NIC) 



Let us now evaluate the EOC on geometric grids which is a particular case of 
non-smooth grids. The numerical computations have been performed on the same 
computational domain as is for the case of geometric grids considered with the normal 
initial condition. The numerical results are presented in Table 15.71 which shows once 
again the convergence of second order. 

The EOC has been computed once more on locally uniform, oscillatory and ran- 
dom meshes for the case mentioned above. The computational domain for locally 
uniform and random meshes is identical as for the previous case. However, we per- 
form the computations on an oscillatory mesh using the same domain as is for the 
uniform mesh. The numerical result are demonstrated in Tables 15. 81 15.91 and 15.101 
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These tables show that we acquire the convergence of first order on locally uniform 
mesh while the fixed pivot technique is zero order convergent on oscillatory and ran- 
dom meshes. 



Table 5.7 
Non-uniform smooth grids (EIC) 



(a) K(x,y) =k (x + y) (b) K(x, y) = koxy 



GP 


Error L\ 


EOC 


Error L\ 


EOC 


60 


7.1E-3 




6.3E-3 




120 


1.8E-3 


1.97 


1.6E-3 


1.98 


240 


4.5E-4 


2.00 


4.0E-4 


2.00 


480 


1.1E-4 


2.03 


1.0E-4 


2.00 



Table 5.9 
Oscillatory grids (EIC) 



(a) K(x,y) = ko{x + y) (b) K(x, y) = koxy 



GP 


Error L\ 


EOC 


Error L\ 


EOC 


60 


0.0554 




0.0295 




120 


0.0532 


0.05 


0.0298 


-0.01 


240 


0.0539 


-0.01 


0.0279 


0.09 


480 


0.0524 


0.04 


0.0256 


0.12 



Table 5.8 
Locally uniform grids (EIC) 



(a) K(x,y) = k (x + y) (b) K(x, y) = koxy 



GP 


Error L\ 


EOC 


Error L\ 


EOC 


60 


0.0303 




0.0145 




120 


0.0156 


0.96 


7.1E-3 


1.04 


240 


7.7E-3 


1.02 


3.3E-3 


1.08 


480 


3.8E-3 


1.03 


1.6E-3 


1.06 



Table 5.10 
Non-uniform random grids (EIC) 



(a) K(x,y) = ko{x + y) (b) K(x, y) = koxy 



GP 


Error L\ 


EOC 


Error L\ 


EOC 


60 


0.0246 




0.0162 




120 


0.0292 


-0.25 


0.0204 


-0.34 


240 


0.0319 


-0.13 


0.0222 


-0.12 


480 


0.0380 


-0.25 


0.0232 


-0.07 



6. Conclusions. A detailed study on the convergence of the fixed pivot tech- 
nique is given for solving aggregation PBE. It is ascertained that the technique is 
second order convergent on uniform and non-uniform smooth meshes. However, it is 
only first order convergent on locally uniform meshes. Finally, the scheme is exam- 
ined closely on oscillatory and non-uniform random meshes and it is observed that the 
technique is not convergent. Furthermore, all observations are validated numerically. 
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